home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / PINVS.PAS < prev    next >
Pascal/Delphi Source File  |  1991-05-01  |  2KB  |  73 lines
  1. PROCEDURE pinvs(ie1,ie2,je1,jsf,jc1,k: integer; VAR c: glcarray;
  2.       nci,ncj,nck: integer; VAR s: glsarray; nsi,nsj: integer);
  3. (* Programs using routine PINVS must define the types
  4. TYPE
  5.    glcarray = ARRAY [1..nci,1..ncj,1..nck] OF real;
  6.    glsarray = ARRAY [1..nsi,1..nsj] OF real;
  7. in the main routine. *)
  8. CONST
  9.    zero=0.0;
  10.    one=1.0;
  11.    nmax=10;
  12. VAR
  13.    js1,jpiv,jp,je2,jcoff,j,irow,ipiv,id,icoff,i: integer;
  14.    pivinv,piv,dum,big: real;
  15.    pscl: ARRAY [1..nmax] OF real;
  16.    indxr: ARRAY [1..nmax] OF integer;
  17. BEGIN
  18.    je2 := je1+ie2-ie1;
  19.    js1 := je2+1;
  20.    FOR i := ie1 TO ie2 DO BEGIN
  21.       big := zero;
  22.       FOR j := je1 TO je2 DO IF (abs(s[i,j]) > big) THEN  big := abs(s[i,j]);
  23.       IF (big = zero) THEN BEGIN
  24.          writeln('pause in routine PINVS');
  25.          writeln('singular matrix - row all 0'); readln
  26.       END;
  27.       pscl[i] := one/big;
  28.       indxr[i] := 0
  29.    END;
  30.    FOR id := ie1 TO ie2 DO BEGIN
  31.       piv := zero;
  32.       FOR i := ie1 TO ie2 DO BEGIN
  33.          IF (indxr[i] = 0)  THEN BEGIN
  34.             big := zero;
  35.             FOR j := je1 TO je2 DO BEGIN
  36.                IF (abs(s[i,j]) > big)  THEN BEGIN
  37.                   jp := j;
  38.                   big := abs(s[i,j])
  39.                END
  40.             END;
  41.             IF (big*pscl[i] > piv)  THEN BEGIN
  42.                ipiv := i;
  43.                jpiv := jp;
  44.                piv := big*pscl[i]
  45.             END
  46.          END
  47.       END;
  48.       IF (s[ipiv,jpiv] = zero) THEN BEGIN
  49.          writeln('pause in routine PINVS');
  50.          writeln('singular matrix'); readln
  51.       END;
  52.       indxr[ipiv] := jpiv;
  53.       pivinv := one/s[ipiv,jpiv];
  54.       FOR j := je1 TO jsf DO s[ipiv,j] := s[ipiv,j]*pivinv;
  55.       s[ipiv,jpiv] := one;
  56.       FOR i := ie1 TO ie2 DO BEGIN
  57.          IF (indxr[i] <> jpiv)  THEN BEGIN
  58.             IF (s[i,jpiv] <> zero)  THEN BEGIN
  59.                dum := s[i,jpiv];
  60.                FOR j := je1 TO jsf DO s[i,j] := s[i,j]-dum*s[ipiv,j];
  61.                s[i,jpiv] := zero
  62.             END
  63.          END
  64.       END
  65.    END;
  66.    jcoff := jc1-js1;
  67.    icoff := ie1-je1;
  68.    FOR i := ie1 TO ie2 DO BEGIN
  69.       irow := indxr[i]+icoff;
  70.       FOR j := js1 TO jsf DO c[irow,j+jcoff,k] := s[i,j]
  71.    END
  72. END;
  73.